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Abstract 

The impressive progress of the kinetic schemes in the solution of gas dynamics problems and the development of 
effective parallel algorithms for modern high performance parallel computing systems led to the development of 
advanced methods for the solution of the magnetohydrodynamics problem in the important area of plasma physics. 
The novel feature of the method is the formulation of the complex Boltzmann-like distribution function of kinetic 
method with the implementation of electromagnetic interaction terms. The numerical method is based on the explicit 
schemes. Due to logical simplicity and its efficiency, the algorithm is easily adapted to modern high performance 
parallel computer systems including hybrid computing systems with graphic processors. 
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1. Introduction 

The tremendous progress in the development of high performance computing systems, especially expecting drasti- 
cally new exascale computing systems, including the challenges in architecture, scale, power and reliability, gives new 
opportunities for the mathematical modeling of important physical phenomena in the present and future. Nevertheless 
the complexity of the challenges in science and engineering continues to outpace our ability to adequately address 
them through impressively growing computational power 

A feature of the present is that the development of technologies and computer systems architecture are well ahead 
of software development. The software problems are primarily associated with the complexity of the algorithms 
adaptation for the differential equations of mathematical physics to high performance computing systems architecture. 
In particular they refer to one of the important requirements as the accuracy in combination with the correctness of 
the initial mathematical models. Another requirement for the methods is their logical simplicity and high efficiency at 
the same time. The numerical algorithms should be simple and transparent from a logical point of view. 

One of the important directions to overcome these problems is the development of a nontraditional approach to 
initial mathematical models and computational algorithms. In the present study for the solution of the multidimen- 
sional gas dynamics and magnetohydrodynamics problems kinetic difference scheme is proposed. It is convenient 
from the physics point of view to define the gas dynamics and magnetohydrodynamics quantities from close relations 
between the kinetic and gas dynamics description of physics processes |[Tl|2l- 

Another aspect is the study of the explicit finite difference schemes, which seem to be preferable for future high 
performance parallel computing, especially in terms of their simplicity and well adaptability to parallel program 
realization, including hybrid high performance parallel computing systems. The weakness of explicit schemes is 
a strictly limited time step that ensures computational stability. This restriction becomes critical with the growing 
number of nodes and the reduction in the step of a spatial mesh. The advanced explicit kinetic finite difference 
schemes have a soft stability condition giving the opportunity to enhance the stability and to use very fine meshes 131 . 
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The mentioned aspects are used for the development of the framework for the study of the dynamics of the 
conducting gas media in strong magnetic fields at high performance parallel computing systems. 



2. Theoretical Issues 

2.1. Gas Dynamics Processes 

The kinetic theory describes the gas dynamics by the Boltzmann differential equation through the evolution of the 
distribution function / (x, ^,t) iSJ: 

^^y^ +^-V/(x,g.O^C(/) (1) 
at 

where C if) is a nonlinear integral operator which describes the collisions between gas molecules. 

This evolution equation follows naturally from the relations between the kinetic and the gas dynamics description 
of continuous media. The macroscopic observables such as density, momentum, energy flux as a function of x and t 
are obtained from the moments of the distribution function with respect to the macroscopic velocity. The evolution 
equations for these gas dynamics quantities are obtained by integrating Eq. ^ over molecular velocities ^ with 
summational invariants m,m^,jm^^. The computational interest in kinetic formulations of the gas dynamics is high 
due to the linearity of the differential operator on the left side of Eq. ([T]i. Nonlinearity is confined by the collision 
term, which is generally local in x and t. 

An important feature is that the collision integral vanishes in the equilibrium state when the local Boltzmann 
distribution function / is a Maxwellian: 

This leads to the use of this model for numerical methods and possible generalizations in order provide a natural 
kinetic description of the system of conservation laws. This approximation is sufficient for the gas dynamics processes 
and is called the kinetic approach |[T1 . 

2.2. Electromagnetic Processes 

In Is) it was shown that electromagnetic fields do not destroy the validity of the Boltzmann equation and this 
opened the way for the implementation of the electromagnetic term in the Boltzmann-like distribution function. From 
the vector nature of the electromagnetic interaction, the distribution function should taking to account the vector 
behavior and provide correct formulation for the evolution of the magnetic field, i. e. the magnetic field should be 
generally defined as the momentum of the Boltzmann-like distribution function. 

A few useful attempts to formulate the vector Boltzmann-like distribution function can be found in |l7][8l|9], but 
physical meaning was not clear defined. 

We propose an evaluation of the electromagnetic processes in the context of the distribution function, taking to 
account the axial nature of the magnetic field. The electromagnetic field is considered as a complex vector field as 
proposed in ||6l: 

F = E + (3) 

For the purposes of magnetohydrodynamics, the efifect, which a magnetic field exerts on a certain volume, is 
obtained by integrating the electromagnetic stress tensor over the surface of that volume and the correspondent prop- 
agation velocity can be defined as a complex vector of velocity: 

^ em ~ ^em em (4) 

At first approximation the term defined of electric forces could be neglected and the magnetic term could be 
defined through the tension of the magnetic field line and shows a similarity to the Alfven wave mechanism: 

y^em = —p (5) 
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2.3. Proposed Distribution Function for MHD 

Using the above definitions we define the local complex Boltzmann Maxwellian distribution function of magne- 
tohydrodynamics with drift velocity u in magnetic field B at the equilibrium: 

The first term on the right-hand side of (j6]l includes the internal energy and the second term is the magnetic 
field energy. The hydrodynamics observables are real scalars and vectors. The complex components include the 
dynamics of the macroscopic observables introduced by the evolution of the magnetic field, keeping their specific 
pseudo-vectorial nature. 

The magnetogasdynamics observables are obtained as integrals of the distribution function (|6| with the sum- 
mational invariants (m,m^, jm^^,m^*). The integration is performed on the path y with respect to the molecular 
velocities ^ in the complex plane correspondent to the value of B in the imaginary space. The relations are obtained 
for the real and imaginary terms; 

p(x,f)= r mfMd^^ (7a) 

Jy 

u(x,r) = — ^ r m^f^d'^ (7b) 

P (X, t) Jy 

E{x,t)^ j^m^^fMd'^ (7c) 

B(x,f) = r mefMd^ (7d) 

^jp (x, f) Jy 

The proposed complex Boltzmann Maxwell like distribution function contains the hydrodynamics terms and the 
electromagnetic terms. Thus by using this distribution function to calculate the mass, momentum, energy and magnetic 
field fluxes, most of the electromagnetic contributions are calculated directly, i.e. one does not have to solve the 
hydrodynamics and magnetic force components separately or differently, as will be shown below. 



3. Ideal MHD System of Equation 



To provide the first step of the formulation of the MHD conservation laws equation, the equilibrium state is 
considered with the proposed distribution function. The MHD system of equations is obtained by the integration of 
([T]l with vanishing collision integral with the summational invariants following the definition in (j7]i: 



ri .df ri 
li'^^Tt^Xi' 

{ mC^ + — [ mf^ ■ Vfd^^ = 
yfp Jy at ^ Jy 



,nf^-Vfd^^O 



The result obtained, set of Eq. (j9]), is the ideal magnetohydrodynamics system of equations: 
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In addition an equation for V ■ B is obtained as the imaginary part of the path integral of the summational invariant 
(m) with respect to the velocities ^: 



Jy ot ^ Jy 



(10) 



p-0 (11) 

oxi 



4. Kinetic MHD Finite Difference Sclieme 



The model of the kinetic differential schemes is based on the discrete model of evolution of the distribution 
function. Kinetic schemes are obtained directly from the Boltzmann kinetic equation by using the principle of total 
approximation. 

Consider the local volume of gas (cell /) with the distribution function in time By using the splitting method of 
particle flow for the cell /, the evolution of the distribution function by first order differential scheme for the kinetic 
Boltzmann equation can be written as: 

- f _ JL -fU \^\ fL-v/+fL ,,,, 

At ~^ 2Ax 2 Ax ' ^ ^ 

As mentioned before the collision of particles leads to the establishment of the equilibrium state which is ade- 
quately described by the single-particle Maxwell distribution function with vanishing of the collision integral in the 
right part of the balance relations. The time evolution of the distribution function can be represented as the time 
evolution of the local Maxwellian distribution function in discrete moments: 

• at time f on each cell, the locally constant one-particle Maxwellian distribution function is defined: 

(13) 



pm '^^ m 



B ' 



where the magnetohydrodynamics parameters p,u,r,B are not varied on the cell. 

• during the time interval At = f^^' - t-' coUisionless processes of the gas dynamics occurs, 

• at time f^^' the distribution function is instantaneously maxwellised 

• for the time f-'^^ these processes are repeated. 

The kinetic difference scheme in this case can be written: 

■fi ~ fi,M _ >.fLl~ fi-\ ^x\^\ f i+\,M ^ i,M f i-\,M .-..^ 

Af ~^ 2Ax 2 'Ax^ 

or in more general form for the multidimensional case: 



where: 

Act, is the surface element Ax^ Ax,,, perpendicular to the direction x,, 

is the value of the distribution function at the surface cr between the two volume elements /, and 



dx, 



the distribution function derivative at the surface cr between the two volume elements. 



The sum in Eq. ( 15 1 is extended to the 6 surface elements at the boundary of the 3-dimensional rectangular volume 
element. 



cesses can be obtained by integrating the balance relation (15 1 with the summational invariants m,ni^, ^in^^,ni^*. 



The kinetic scheme of the conservation laws of the macroscopic observables for 3D magnetohydrodynamics pro 
ises can be obtained by integrating the ba 
using the same integration rules as in Eq. (j9]i: 

Ax; 
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where p - 



2p + B^ 



,i,k^l...3 



In addition to Eq. ( 16 1 the condition is obtained as the complex part of the path integral of the summational 



invariants (m) with respect to the molecular velocities ^ in the complex plane: 

(Bih, = ^ [Bih,,, 



(17) 



Dissipative terms appear in the time evolution of the magnetic field which does not preserve the condition V ■ B = 
and require a specific treatment. 



5. Kinetic Quasi MHD Equations 

The kinetic quasi magnetohydrodynamics system of equations is closely related to the kinetic scheme and repre- 
sents a differential form notation for the numerical algorithms. 
The balance relation in Eq. ( [T5| ) can be rewritten as: 



and using the Gauss-Ostrogradsky formula it is possible to transform Eq. ( 18 1 to the differential form 

dt y^-^^i 2 dxi dxk 



'^■{uL)-~,^^>^^fL (19) 



Here the quasi magnetohydrodynamics system of equation involves explicitly two t parameters. Hydrodynam- 
ics processes are introduced by the quantity t that corresponds to the time of free distance flight of particles, or the 
characteristic time of particle collisions. By analogy the quantity t,„ is introduced as the characteristic time of propa- 
gation of magnetohydrodynamics by electromagnetic processes. The characteristic time values r and t,„ are defined 
respectively for hydrodynamics and electromagnetic processes: 

Ax,- Ax/ 
r = a— r„, = a,„— (20) 

where: 

Ax„ is the size of the computational cell, 

Ch, Cm are the sound speed and Alphen speed in the computational cell. 

The introduction of the physical meaning of characteristic times t and t,„ provides an important contribution to 
the understanding of the processes and the simplification of the numerical scheme. 

The evolution equations for the gas dynamics parameters and for the magnetic field are obtained from Eq. ( 19 1 by 
integration with the summation invariants (pi,^) - m,m^, ^m^^,m^* over the molecular velocities, under the assump- 
tion: 

j f^'^i^)d€^ j f^'cl>{^)d^ (21) 

The integration is performed as in Eq. (|9]l and Eq. ( 16 1. The gas dynamics and magnetic field quantities are obtained 
respectively as the real and imaginary path of the integral in the complex space. 

The compact form of the kinetic quasi magnetohydrodynamics system of equations can be written as: 

p^+' 'pJ d dwi 

— + '^P^i ^ ^ (22a) 

Af OX; OX/ 

- — ' , ' + ^n,i = —nl + —WiUk (22b) 
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^—^ + ^Mf, = ^nf (22d) 
At dxk dxk 



The left-hand part of the system of Eq. ( p2| corresponds to the Euler system of equations. The right-hand of the 
kinetic quasi MHD Eq. ( [22] i contains dissipative terms. In comparison with other methods, the dissipative terms are 
obtained not by phenomenology with some assumption about magnetohydrodynamics processes but in consistency 
with the diff'erence scheme of the Boltzmann equation. 

n,i is the momentum flux density tensor for a perfect gas in magnetic field: 

/ B^\ 

n,A- ^\p + — \5ik+ puiUk - BiBk (23) 



Fi is the heat transfer flux of a perfect gas in magnetic field: 



■ p + —\ Ui - BiUkBk 



(24) 



Mfi^ is the asymmetric product between velocity u and magnetic field flux B: 



Ml = UkBi - UiBk (25) 



The right hand part of the system of Eq. ( 22 1 includes the dissipative terms: 



6 



T_d_ 

2 dxk 



2 dxk 



(26a) 



dui duk 



T 

^2 



B- 



2 du„, 

3 ox,„ 



dui 



B' 



\ duk I B- 



— 5mk - B,„Bk I H I —5im - BiB„,j — | —6i, 



dx„ 



dui diik dii„ 
Bm I -Ba J, Bi h B„- — 6ik 

OXn! OXi^, OXm 



piliU,, 



2 

duk 
dx„, 



dx,„ 



dp d B- d 
+ Ui— — Ui- — B„,Bi, 
dxk oxk 2 dx„, 



dp du„ 



Sik 



du,n dUk 
-BiBk- — + BiB„,- BiU 



2 du,„ 
dXf^ 
dBk 



dUfi dBfj 

' Bn Bfii — + Bf^ Ufji — 

OXm OXm 



dx„. 



dx„, 



dxn. 



Sik 



-BkBi 



du„ 

dx„. 



+ BkBn 



dui 

' dx,„ 



- BkU 



dBi 

Oxff, 



(26b) 



5 d p 
2 ^ dxi p 



2 



B^ 



2 

p5ik 



Sik - DjDk 

B^ 



j dxkp 



pUjUk 



2 

d 3p 
dxk 2 p 



5ik - BiB, 



V 



-pUiUk 

P + 
P + 



d_lf_ 
j dxk 2p 

pUiUk [p + B-) 

B- d I 1 d B^ 
2 dxk p p dxk 2 

B^MdBi^ _ dBk \ 

dxk dxi j 

d I d 1 
Bk 



P + 



B^\ d BiBk BjBk d B^ 



d 1 



dxk p 

— M,D — 

dxk p dxk 
du,„ 



BiBn 



-Uk 



dxk 



B, 



dxk p 
1 

- -BkB„ 
P 



dxi p 
dBi 



UkB,, 



du, 

dx,n 



p dxk 2 

I- z UiBm I B„t 
2 [ \ dxk 

1 dp Id B- 

p dxi„ p dx„, 2 



duk 
' dx,„ 



duk du,„ dB,„ 

t>k^ — + Uk— — 

dxk dxk 



l_d_ 
pdx, 



B,„Bk 



1 dBk 

-BiBm- 

p dx,„ 

^ duk Bi dp Bi d B^ 

BiUm ^ + — + — -7. 

dx„, p dxk p dxk 2 



dum dui dBi 

UkBi- UkB„,- + UkUm- — 

OXfyi OXm OXm 



Bi d 
P dx„ 



du„ duk 

-UitSk- + UitS„,- 

oXffi oXffi 



dBk 

' dx,n 



BkB„ 



-BkU, 



(26c) 



(26d) 



du^ _ Bkdp _ Ihd_lf_ ^ Ih^gg 
' dx„ p dxi p dxi 2 p dx,„ ' " 



The dissipative terms appear because the construction of the quasi magnetohydrodynamics system is based on 
the assumption that the distribution function sHghtly changes over the distance between neighborhood cells, what is 
related to the characteristic times t and t„,. It was shown in 12] that the dissipative terms of the quasi gas dynamics 
system are small in comparison with the convective terms with the condition of cell size equivalent to the free path 
they converge to the viscous terms of the corresponding Navier-Stokes equations. The corresponding dissipative terms 
are associated with real physics processes. An important remark is that in this case the gas dynamics parameters such 
as viscosity and heat conductivity are obtained from the kinetic theory. 



The Navier-Stokes viscosity is identified as the first term of Eq. (26 5): 
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(27) 



where the bulk viscosity component is neglected and the shear viscosity n is related to the gas pressure p and the 
characteristic time t as = ^p. 



The Navier-stokes thermal flux vector is identified as the first term of Eq. (26 1): 



5 d p 
2 axi p 



dxi 



(28) 



with T gas temperature and k thermal coeflicient expressed as A: = ^ f^jP- with Pr Prandl number. 

A similar analysis of the dissipative terms of the electromagnetic processes gives the estimation of their smallness. 
With correct conditions for the size of cells the equation converges to the correct representation of magnetic viscosity. 



The gas resistivity is identified as the first term of Eq. (26 ;) and also appeares as a result of the kinetic formulation; 



p + 



dBi 
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(29) 



with the resistivity 77 = ^(z^ + t) 



6. Computational Algorithm 

The computational algorithm uses a Cartensian, staggered, divergence free mesh configuration. A detailed de- 
scription is presented in ifTOlfTTI . in order to preserve the condition V ■ B = 0. 

Fig. [T] shows the four neighbour to the cell {i,k, I) used in evaluations of the hydrodynamics and electromagnetic 
variables. 

The hydrodynamics observables - mass density, momentum and energy density are defined at the cell center. The 
components of the magnetic field are defined at the face centers of the cells. A duality is established between the 
electric field and the fluxes. This duality is utilized to obtain the electric field at the edges of the computational cell 
through a reconstruction process that is applied directly to the properly upwinded fluxes. The electric field is then 
utilized to make an update of the magnetic fields that preserves the solenoidal nature of the magnetic field and ensures 
that the magnetic field in a magnetohydrodynamics model remains strictly solenoidal up to discretization errors. 



p,prf,E 




Ey 



Ez 
Ey 



■Ex^ I Ex Ez^ 



Figure 1 : 3D Computational Domain 

Generally the explicit numerical scheme is used model, considering that it is perspective for the modern high 
performance computing systems due to the logical simplicity and efficiency of the algorithms. 
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The finite volume method is used to update the conserved observables, mass, momentum and energy, by calculat- 
ing the fluxes of this observables across the cell face. Updating the magnetic field is a more complicated procedure 
and is performed via electric field integration along the edge of the cells, as showh on Fig.[T] The distribution function 
method is proposed in the calculation as described above. 

The explicit scheme is used in the time evolution for the integration of the quasi magnetohydrodynamics system 
of equations. The code uses a variable time step. The time step in an explicit scheme is controlled by a Courant type 
conditions on the time step estimation 121 . 



7. Results of Numerical Modeling 

The computational framework is created on the basis of Fortran 90 and C++, with parallel implementation on MPI. 
The demonstration of the performance of the method is performed on the basis of the solution of the spherical 
expansion problem of ionised gas and the solution of the expansion of an ionised gas in strong magnetic field. 
The simulations are performed for a Cartesian rectangular mesh 100 x 100 x 100 in the physics domain [0,1] . 
The initial conditions consist of a sphere with radius 0. 1 placed in the center of the physical region with pressure 
of 100 in comparison to the overall represented area with pressure 1. For the study of ionized gas in a strong magnetic 
field the uniform magnetic field aligned with the z coordinate is added to the initial conditions. 

Fig. 2]3 present the state of the 3D simulation of the processes for relative time 0.03. On the 3D pictures the 



arrows represent the velocities of the ionised gas and the color represents the density of gas. 3D figures clearly show 
the confinement of the ionized gas in the cylindrical area along z due to the magnetic field. 

Fig. |4|5 1 represent the 2D projections of the density, pressure and kinetic energy of the gas expansion without 
magnetic field and the ID density profile for these condition. 

Fig. 6j7 represents the 2D projections of the density, pressure, magnetic pressure and kinetic energy for the gas 



expansion problem of the ionized gas with initial magnetic field of 5/ y/n and Fig.jsjshows the ID profile of the density 
for these conditions. 

Fig. 9|10 represent the 2D projections of the density, pressure, magnetic pressure and kinetic energy for the gas 



expansion of the ionized gas in the strong magnetic field of 50/ yfn at time 0.01 and Fig. 1 1 represents the ID profile 
of the density for these conditions. 

Similar studies of the problem of spherical explosion including the conditions with magnetic field are presented 
in lfT2ll . The comparison of the results shows a reasonable agreement and will be analysed further. 



8. Conclusions 

A new 3D kinetic algorithm has been developed for the solution of the magnetohydrodynamics problems. The 
novel feature of the method is that the local complex Boltzmann-like distribution function incorporated most of the 
electromagnetic processes terms. The fluxes of mass, momentum and energy across the cell interface as well as the 
magnetic field are calculated by integrating a local complex Boltzmann-like distribution function over the velocity 
space. Thus by using this distribution function to calculate the mass, momentum and energy fluxes, most of the 
electromagnetic contributions are calculated directly, i.e. one does not have to solve the hydrodynamics and magnetic 
force components separately or differently. 

A staggered, divergence free mesh configuration is used for the evaluation of the electromagnetic behaviour. 

Numerical examples demonstrate that the proposed method can achieve high numerical accuracy and resolve 
strong shock waves of the magnetohydrodynamics problems. 
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Figure 2: 3D view of tlie gas expansion witiiout magnetic field 
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Figure 9: 2D gas density and 2D gas pressure projections in tlie magnetic field 50/ yfn 
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